In this case, we set K_i = 1 and G_j = K ω_j. The mean-field theory is useful to investigate the stationary solutions of the system in the long time limit, and the basic idea is to conduct self-consistency equations and mean-field frequency. In addition to the traditional order parameter Z(t), a weighted order parameter is necessary, which is defined as 4 With this definition, dynamical equation (3) can be rewritten in the mean-field form 5 Here K|D| can be interpreted as an effective coupling strength. Equation (5) reflects that the interactions between all oscillators are proportional to the coupling with coefficient D. To study the steady states of the system, the self-consistent method turns out to be effective. We assume equation (5) approaches to a stationary state in a long time limit. At this point, the collective amplitude |D| is time independent and the mean-field phase ϕ rotates uniformly with a frequency Ω. We introduce the phase difference 6 The equation of motion of the phase difference oscillator can be transformed into 7 In the thermodynamical limit N → ∞, we introduce the probability density distribution function as ρ(φ, ω,t), which denotes the probability distribution of state variables (φ,ω,t) and satisfies the normalization condition and periodic condition ρ(φ + 2π, ω,t) = ρ(φ,ω,t). Consequently, equation (7) is equivalent to the following continuity equation: 8 The stationary solutions of Eq. (8) can be divided into phase-locked case (|ω − Ω| < K|D|) and drifting case (|ω − Ω| > K|D|). Considering both contributions of the phase-locked and drifting oscillators to D, and separating the real and imaginary parts of the amplitude of D, we obtain the following self-consistency equations corresponding to the real and imaginary parts of |D|: 9 10 Here Θ(x) is the Heaviside function that makes it positive in the radical formula. Theoretically, all steady states of the system can be calculated through Eqs. (9) and (10). In particular, in the limit case (|D| → 0, Ω → Ω_c), the critical point K_c for phase tradition can be obtained analytically as 11 and the critical mean-field frequency Ω_c satisfies the balance equation 12 The symbol P means the principal-value integral along the real ω.^[19]

Another way to get the critical points K_c is to perform the linear stability analysis of the incoherent state. With the density distribution function F(θ,ω,t) = ρ(φ + Ω t, ω,t), the velocity in the continuity equation should be written as 13 The order parameter D(t) can be expressed as 14 Let 15 be the n-th Fourier coefficient of F(θ,ω,t). Substituting it into the continuity equation, we can obtain 16 which must be solved self-consistently with the equation, 17 It is obvious that the incoherent state (|D| = 0) implies that Z_n = 0 for all n (n ≠ 0). Linearizing Eq. (16) around this fixed point, we have 18 The linearized equation (18) can be transformed into eigenvalue equation 19 It is convenient to investigate eigenvalue equation into Cartesian coordinate system by setting λ = x + iy. The sign of x determines the stability of the incoherent state. As a well known result, the eigenvalue λ does not exist when the coupling strength K is small enough provided that ω g(ω) has no singularity. Once K > K_c, λ appears with x > 0 and the incoherent state loses its stability. The critical point K_c can be obtained by imposing the conditions x → 0⁺, y → Ω_c, which has the same formula as Eqs. (11) and (12).

In this situation, we set K_i = Kω_i and G_j = 1 that corresponds to the conformist and contrarian model studied in Ref. [34]. Conduct self-consistency equations for r and mean-field frequency Ω, which are 20 21 where sgn(ω) is the sign function. In the limit case r → 0⁺, Ω → Ω_c, considering the Taylor series of Eqs. (20) and (21) and preserving the first order term, one obtains the same results as Eqs. (11) and (12).

Following the same procedure, the linear stability analysis of the incoherent state can be performed by changing the variable D_n = ω Z_n. Once again, the equation of eigenvalue λ = x + iy for δ Z₁ has the same result as Eq. (19). Letting x → 0⁺, y → Ω_c, we can obtain the same formula for K_c as Eqs. (11) and (12). Therefore, we conclude that the position of the heterogeneity of the coupling strength (in or out) does not influence the critical point for the system where a phase transition takes place and the incoherent state loses its stability in the limit N → ∞.

It is clear that the linearized equation can be solved by making Laplace transform δ Z₁(t,ω) → δ Z₁(s,ω), we have 22 Solving it, one obtains 23 Integrating it with respective to ω for both sides, we obtain 24 Considering the definition of δ Z(t), equation (24) becomes 25 After simplification, we have 26 For convenience, we assume that δ Z₁(t = 0,ω) is a constant, i.e., δ Z₁(t = 0, ω) = 1, which means that the initial phases of all oscillators are identical, and g(ω) follows the Lorentzian distribution 27 where γ denotes the width of the distribution, and Δ stands for the center. We have 28 29 The perturbed order parameter δZ(t) can be determined through the inverse Laplace transform of δZ(s), which reads 30 Obviously, the critical coupling strength is 31 which is consistent with the result obtained by the mean-field theory and linear stability analysis as shown in Fig. 1. Notice that the order parameter decays in the exponential form while the factor (Δ K/2 − γ, Δ + K γ/2) corresponds to the resonance pole in Ref. [19]. Moreover, it is consistent with the result in Ref. [34].

Fig. 1

	Figure Option ViewDownloadNew Window
	Fig. 1 The critical point for phase transition Kc vs. parameters in heterogeneous case. The natural frequency obeys (a), (b) ; (c) , where p represents the height of peak at Δ, with 0 < p < 0.5; (d) g(ω) = 1/2, where ω ∈ (ε −1, ε + 1) and ε > 0 is a bias of uniform distribution.

We emphasize that the mean value Δ controls the property of the system. When Δ > 0, the attractive coupling dominates the system and the incoherent state becomes unstable until K > K_c. While Δ < 0, the repulsive term governs the system, the incoherent state is always stable once K > 0. Particularly, Δ = 0 is a critical case which leads to K_c → ∞.

In this situation, there are two order parameters that should be determined respectively. Considering the variable substitution δ D₁(t = 0,ω) = ω δ Z₁(t = 0,ω), the expression of δ D(s) can be obtained directly from Eq. (26) 32 Applying Laplace transform for Eq. (18), as δ Z₁(t,ω) → δZ₁(s,ω), we can obtain 33 hence 34 Using the same method, we can obtain 35 Considering the expression of δ D(s), we can obtain the final expression of δ Z(s). Assuming δ D₁(t = 0,ω) = a, δ Z₁(t = 0,ω) = b, where a and b are constants (when the initial phases are identical, b = 1 and a is the mean value of the natural frequency), and ω obeys the Lorentzian distribution as Eq. (27). Taking Eqs. (28) and (29) into consideration and making the inverse Laplace transform, we obtain the weighted order parameter δ D(t) as 36 and the original order parameter δ Z(t) as 37

Obviously, we obtain the critical coupling strength as Eq. (31). The above results show that the Landau damping effect is a common phenomenon, as shown in Fig. 2. The real parts of resonance poles control the exponential relaxation rates of the order parameters. We clarify that the weighted order parameter decays in the exponential form that is the same as the out-weighted case while the original order parameter is controlled by δ D(t). It should be pointed out that, the exponential decay of the order parameter is just a special case which is caused by the regularity of g(ω) and identical initial condition. The general decay of the order parameter is possible according to the specifical perturbations and g(ω).^[20–22]

Fig. 2

Figure Option ViewDownloadNew Window

Fig. 2 The decay of δ D(t) and δ Z(t) in the in-frequency-weighted case with different K (K < Kc). The solid lines are the direct numerical solutions of Eq. (3) (here Ki = 1 and Gj = K ωj) and the dashed lines are the fitted curves with exponent as Eqs. (36) and (37). The natural frequency obeys . (a) and (b) K = 1.6, Δ = 0.5, γ = 0.5. (c) and (d) K = 1.4, Δ = 0.5, γ = 0.5. In the numerical simulations, we set the total number of oscillators to N = 100000, and the initial phases for all oscillators are identical.